%% Function to compute the equilibrium in urban areas
function f = UrbanEquilibriumFunction(tau_F_income,pp,pop_C,pop_F,xiH_F,xiH_C,XS,XM,varphi,QA_eq, alphaA,alphaM,alphaS,cAbar,cSbar, Xn,Xr,XrXr,XnXn,f_XrXn_C,f_XrXn_F, L_r_M_only,L_n_M_only,Q_M_only,pSM_M_only);

%% Unpackaging prices   
    p_SU    = pp(1);
    p_M     = pp(2);
   
%% Constructing measures of skills distributions in Urban Regions (U= C + F)    
   f_XrXn_U   =  f_XrXn_C +f_XrXn_F;
   
    Xn_f_XrXn_C   =  XnXn.*f_XrXn_C;
    Xn_f_XrXn_F   =  XnXn.*f_XrXn_F;
    Xn_f_XrXn_U   =  XnXn.*f_XrXn_U;

    Xr_f_XrXn_C   =  XrXr.*f_XrXn_C;
    Xr_f_XrXn_F   =  XrXr.*f_XrXn_F;
    Xr_f_XrXn_U   =  XrXr.*f_XrXn_U;

    factor_wn  =  varphi^( varphi/(1-varphi) )*(1-varphi)*(XM/XS^varphi)^(1/(1-varphi));

    if p_SU/p_M < pSM_M_only,
        wrU     =       varphi*p_M*XM*(L_n_M_only/L_r_M_only)^(1-varphi);
        wnU     =   (1-varphi)*p_M*XM*(L_n_M_only/L_r_M_only)^(-varphi);
    else
        wrU     =    p_SU*XS;
        wnU     =    factor_wn*( p_M /(p_SU)^varphi )^(1/(1-varphi));  % these are the wages implied by the case when both S and M are produced at U
    end

%% Individual Earnings
    yU                              =    max(wrU*XrXr,wnU*XnXn);
    
%% Consumption of Services
    ybar_F                          =    p_SU*cSbar*(1-alphaS)/alphaS + cAbar + p_M*xiH_F;        % income threshold for the consumption of Services, households in slums
    ybar_C                          =    p_SU*cSbar*(1-alphaS)/alphaS + cAbar + p_M*xiH_C;        % income threshold for the consumption of Services, households in cities
    ICS_F                           =    zeros(size(XrXr));  % Indicator, positive consumption, households in slums
    ICS_C                           =    zeros(size(XrXr));  % Indicator, positive consumption, households in cities
    ICS_F(yU>=ybar_F)               =    1;
    ICS_C(yU>=ybar_C)               =    1;
% Aggregation: Consumption of Services
    yTop_F                  =    trapz( Xr,trapz(Xn, yU.*f_XrXn_F.*ICS_F));    % Slums
    fTop_F                  =    trapz( Xr,trapz(Xn,     f_XrXn_F.*ICS_F));
    yTop_C                  =    trapz( Xr,trapz(Xn, yU.*f_XrXn_C.*ICS_C));    % Cities
    fTop_C                  =    trapz( Xr,trapz(Xn,     f_XrXn_C.*ICS_C));

   
% Aggregation: Supplies of Skills and Employment    
    IOCC_r_U                        =    zeros(size(XrXr));  % indicator, occupation choices
    IOCC_r_U(wnU*XnXn < wrU*XrXr)   =    1;
    LrU                             =    trapz( Xr,trapz(Xn, Xr_f_XrXn_U .*IOCC_r_U));      %  effective labor in routine tasks, urban areas
    ErU=    trapz( Xr,trapz(Xn,    f_XrXn_U .*IOCC_r_U));      %  employment in routine tasks, urban areas
    LnU=    trapz( Xr,trapz(Xn, Xn_f_XrXn_U .*(1-IOCC_r_U)));  %  effective labor in non-routine tasks, urban areas
    EnU=    trapz( Xr,trapz(Xn,    f_XrXn_U .*(1-IOCC_r_U)));  %  employment in non-routine tasks, urban areas
    if p_SU/p_M<pSM_M_only,
        QM=   Q_M_only;
        QSU=   0; 
    else
        LrM                                       =   (varphi*XM*p_M/XS/p_SU)^(1/(1-varphi))*LnU;
        QM  =   XM*(LrM)^varphi*(LnU)^(1-varphi);
        QSU =   XS*(LrU-LrM);          
    end
 
%% Equilibrium Pricing Conditions
    PP_SU1  =    alphaS*( yTop_F*(1-tau_F_income)+yTop_C - fTop_C*(cAbar+p_M*xiH_C) - fTop_F*( cAbar+p_M*xiH_F ) )/( QSU +(1-alphaS)*cSbar*(fTop_F+fTop_C) );
    PP_MM1  =    alphaM/alphaA*( QA_eq - cAbar )/( QM -xiH_C*pop_C-xiH_F*pop_F );
     
%% Distance (Euclidean)
    f = (PP_SU1-p_SU)^4 + (PP_MM1-p_M)^4;  
 
 
